//
//		Author: Alexandros Sigalas
//
//      This program is free software; you can redistribute it and/or modify
//      it under the terms of the GNU General Public License as published by
//      the Free Software Foundation; either version 2 of the License, or
//      (at your option) any later version.
//      
//      This program is distributed in the hope that it will be useful,
//      but WITHOUT ANY WARRANTY; without even the implied warranty of
//      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//      GNU General Public License for more details.
//      
//      You should have received a copy of the GNU General Public License
//      along with this program; if not, write to the Free Software
//      Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
//      MA 02110-1301, USA.

#include <stdlib.h>
#include <assert.h>
#include <string.h>
#include <math.h>
#include "globals.h"
#include "mathutils.h"
#include "polygon.h"
#include "material.h"


polygon pl_init(const uint count)
{
    polygon pl= (polygon)calloc(1, sizeof(struct Polygon));
    assert(pl);
    pl->points= (struct PL_Vertex*)malloc(count * sizeof(struct PL_Vertex));
    assert(pl->points);
    pl->size= count;
    return pl;
}

polygon pl_copy(const polygon pl, void (*cpy)(vec * dest, const vec * const src), const int copy_holes)
{
    assert(pl);
    polygon cp= (polygon)malloc(sizeof(struct Polygon));
    assert(cp);
    memcpy(cp, pl, sizeof(struct Polygon));
    cp->points= (struct PL_Vertex*)malloc(pl->size * sizeof(struct PL_Vertex));
    assert(cp->points);
    uint i;
    if (cpy){
        for (i= 0; i < pl->last; i++){
            cpy(&cp->points[i].co, &pl->points[i].co);
            cp->points[i].flags= pl->points[i].flags;
        }
        cpy(&cp->bb.min, &pl->bb.min);
        cpy(&cp->bb.max, &pl->bb.max);
    }
    else{
        memcpy(cp->points, pl->points, pl->last * sizeof(vec));
    }
    if (copy_holes)
        cp->holes= pll_copy(pl->holes, cpy);
    return cp;
}

void pl_kill(polygon *pl_pt)
{
	assert(pl_pt);
	polygon pl= *pl_pt;
	if(pl){
		if (pl->points)
			free(pl->points);
		if (pl->holes)
			pll_killall(&pl->holes);
		free(pl);
	}
}

static void pl_check_bounds(polygon pl, const vec * const p)
{
    if (pl->last > 1)
        check_bounds(&pl->bb, p);
    else{
        v_copy(&pl->bb.min, p);
        v_copy(&pl->bb.max, p);
    }
}

static void pl_recalc_bounds(polygon pl)
{
    assert(pl);
    v_copy(&pl->bb.min, &pl->points[0].co);
    v_copy(&pl->bb.max, &pl->points[0].co);
    uint i;
    for (i=1; i < pl->last; i++)
        check_bounds(&pl->bb, &pl->points[i].co);
}

static polygon pl_check_pointer(polygon *pl_pt){
	assert(pl_pt);
	polygon result= *pl_pt;
	if (!result)
		result= *pl_pt= pl_init(10);
	return result;
}

void pl_expand(polygon pl, const uint extra_space)
{
    assert(pl);
    pl->points= (struct PL_Vertex *)realloc(pl->points, pl->size += extra_space);
    assert(pl->points);
}

struct PL_Vertex *pl_append(polygon *pl_pt, const vec * const co, int flags)
{
	assert(pl_pt);
	polygon pl= *pl_pt;
    assert(pl);
    assert(co);
    if (pl->last > pl->size)
        pl_expand(pl, pl->size);
    struct PL_Vertex *vx= pl->points + pl->last++;
    vx->co= *co;
    vx->flags= 0;
    pl_check_bounds(pl, co);
    return vx;
}

struct PL_Vertex *pl_append_point(polygon *pl_pt, const vec * const p)
{
	
	assert(p);
	polygon pl= pl_check_pointer(pl_pt);
	if (!pl)
		pl= *pl_pt= pl_init(10);
    if (pl->last > pl->size)
        pl_expand(pl, pl->size);
    struct PL_Vertex *vx= pl->points + pl->last++;
    vx->co= *p;
    pl_check_bounds(pl, p);
    return vx;
}

struct PL_Vertex *pl_append_vertex(polygon *pl_pt, const struct PL_Vertex * const pt)
{
    assert(pt);
    polygon pl= pl_check_pointer(pl_pt);
    if (pl->last > pl->size)
        pl_expand(pl, pl->size);
    struct PL_Vertex *vx= pl->points + pl->last++;
    *vx= *pt;
    pl_check_bounds(pl, &pt->co);
    return vx;
}

struct PL_Vertex *pl_insert_point(polygon pl, uint pos, const vec * const p)
{
    assert(pl);
    assert(p);
    if (pl->last > pl->size)
        pl_expand(pl, pl->size);
    uint i= pl->last++;
    uint j= pl->last;
    for (; i >= pos; j=i--)
        pl->points[j]= pl->points[i];
    struct PL_Vertex *vx= pl->points + i;
    vx->co= *p;
    vx->flags= 0;
    pl_check_bounds(pl, p);
    return vx;
}

struct PL_Vertex *pl_insert_vertex(polygon pl, uint pos, const struct PL_Vertex * const p)
{
    assert(pl);
    assert(p);
    if (pl->last > pl->size)
        pl_expand(pl, pl->size);
    uint i= pl->last++;
    uint j= pl->last;
    for (; i >= pos; j=i--)
        pl->points[j]= pl->points[i];
    struct PL_Vertex *vx= pl->points + i;
    *vx= *p;
    pl_check_bounds(pl, &p->co);
    return vx;
}

struct PL_Vertex *pl_pop(polygon pl)
{
    assert(pl);
    if (pl->last){
        struct PL_Vertex *result= (struct PL_Vertex*)malloc(sizeof(struct PL_Vertex));
        pl->last--;
        struct PL_Vertex *last= pl->points + pl->last;
        v_copy(&result->co, &last->co);
        result->flags= last->flags;
        return result;
    }
    return NULL;
}

struct PL_Vertex *pl_remove(polygon pl, const uint idx)
{
    assert(pl);
    if (idx < pl->last){
        struct PL_Vertex *result= (struct PL_Vertex*)malloc(sizeof(struct PL_Vertex));
        memcpy(result, pl->points+idx, sizeof(struct PL_Vertex));
        uint i=idx;
        struct PL_Vertex *vx= pl->points + i++;
        for(;i < pl->last; i++)
                *vx++=pl->points[i];
        pl->last--;
        pl_recalc_bounds(pl);
    }
    return NULL;
}

struct PL_Vertex *pl_remove_vertex(polygon pl, const struct PL_Vertex * const vx)
{
    assert(pl);
    assert(vx);
    uint i;
    for (i=0; i < pl->last; i++){
        if (pl->points + i == vx)
            return pl_remove(pl, i);
    }
    return NULL;
}

void pl_extend_verts(polygon pl, const struct PL_Vertex * const points, uint count)
{
    assert(pl);
    assert(points);
    if (pl->last + count > pl->size)
        pl_expand(pl, pl->size + count );
    uint i;
    for (i=0; i < count; i++)
        pl->points[pl->last++]= points[i];
}

void pl_extend_points(polygon pl, const vec * const points, uint count)
{
    assert(pl);
    assert(points);
	struct PL_Vertex *tmp;
    if (pl->last + count > pl->size)
        pl_expand(pl, pl->size + count );
    uint i;
    for (i=0 ; i < count; i++){
		tmp= pl->points + pl->last++;
		tmp->flags= 0;
		tmp->co= points[i];
	}
}

int pl_add_hole(polygon pl, const polygon hole)
{
    //TODO:fix this
    pll_append(&pl->holes, hole);
    return 1;
}

int pl_rel_point_2d(const polygon pl, const vec * const p, enum axis_pair ax, const int do_holes)
{
    assert(pl);
    assert(p);
    if (!point_in_bounds(&pl->bb, p))
        return OUT;
    uint i;
    if (do_holes && pl->holes){
        int rel;
        for (i=0; i < pl->holes->last; i++){
            if ((rel= pl_rel_point_2d(pl->holes->polys[i], p, ax, 0)) != OUT){
                return (rel == IN) ? OUT : ON;
            }
        }
    }
    uint j, wn= 0;
    double x, y, x0, y0, x1, y1, ux, uy;
    int _x= (ax == YZ) ? 1 : 0;
    int _y= (ax == XY) ? 1 : 2;
    double *v0= (double*)p, *v1;
    x= *(v0 + _x);
    y= *(v0 + _y);
    for (i=(pl->last - 1), j=0; j < pl->last; i=j++){
        v0= (double *)&pl->points[i].co;
        v1= (double *)&pl->points[j].co;
        x0= *(v0 + _x);
        y0= *(v0 + _y);
        x1= *(v1 + _x);
        y1= *(v1 + _y);
        ux= x1 - x0;
        uy= y1 - y0;
        if (fabs((x - x0) * uy - (y - y0) * ux) < E){
            if (fabs(ux) > fabs(uy)){
                if ((x0 < x && x < x1) || (x1 < x && x < x0)){
                    return ON;
                }
            }
            else if ((y0 < y && y < y1) || (y1 < y && y < y0)){
                    return ON;
            }
        }

        if (y0 <= y){
            if ((y1 > y) && ((x - x0) * uy - (y - y0) * ux > 0.0)){
                wn += 1;
            }
        }else if((y1 <= y) && ((x - x0) * uy - (y - y0) * ux < 0.0)){
            wn -= 1;
        }
    }
    return (wn) ? IN : OUT;
}

double pl_signed_area_2d(const polygon pl, enum axis_pair ax)
{
    assert(pl);
    double area= 0.0;
    if (pl->last < 3)
        return area;
    uint i, j;
    for (i=(pl->last - 1), j=0; j < pl->last; i=j++)
        area+= v_perp(&pl->points[i].co, &pl->points[j].co, ax);
    area/= 2.0;
    return area;
}

void pl_area_vec(const polygon pl, vec * dest)
{
    assert(pl);
    v_zero(dest);
    if (pl->last < 3)
        return;
    uint i, j;
    vec tmp;
    for (i=pl->last - 1, j=0; j < pl->last; i=j++){
		v_cross(&tmp, &pl->points[i].co, &pl->points[j].co);
		v_add(dest, dest, &tmp);
	}
}

double pl_area_3d(const polygon pl, const int do_holes)
{
    assert(pl);
    double area= 0.0;
    if (pl->last < 3)
        return area;
    vec tmp;
    pl_area_vec(pl, &tmp);
    area= v_length(&tmp);
    if (do_holes && pl->holes){
        uint i;
        for (i=0; i < pl->holes->last; i++){
            pl_area_vec(pl->holes->polys[i], &tmp);
            area-= v_length(&tmp);
        }
    }
    return area;
}

double pl_area_2d(const polygon pl, enum axis_pair ax, const int do_holes)
{
    assert(pl);
    double area= 0.0;
    if (pl->last < 3)
        return area;
    area= fabs(pl_signed_area_2d(pl, ax));
    if (do_holes && pl->holes){
        uint i;
        for (i=0; i < pl->holes->last; i++)
            area-= fabs(pl_signed_area_2d(pl->holes->polys[i], ax));
    }
    return area;
}

enum PL_Order pl_order_2d(const polygon pl, enum axis_pair ax)
{
    if (pl_signed_area_2d(pl, ax) > 0.0)
        return CW;
    else
        return CCW;
}

enum PL_Order pl_order(const polygon pl)
{
    vec n;
    pl_normal(pl, &n);
    if (fabs(n.z) < ZEROLENGTH){
        if (fabs(n.x) > fabs(n.y))
            return pl_order_2d(pl, XZ);
        else
            return pl_order_2d(pl, YZ);
    }
    else
        return pl_order_2d(pl, XY);
}
    
void pl_normal(const polygon pl, vec *normal)
{
    assert(pl);
    if (pl->last < 3)
        v_zero(normal);
    else{
		pl_area_vec(pl, normal);
		v_normal(normal, normal);
	}
}

int pl_rel_plane(const polygon pl, const vec * const pp, const vec * const pn, int* rmap)
{
    int R= 0;
    uint i;
    if (rmap){
        for (i=0; i < pl->last; i++)
            R|= *rmap++= rel_point_plane(&pl->points[i].co, pp, pn);
        
    }
    else{
        for (i=0; i < pl->last; i++)
            R|= rel_point_plane(&pl->points[i].co, pp, pn);
    }
    return R;
}

int pl_self_intersects(const polygon pl, enum axis_pair ax)
{
    vec diffs[pl->last];
    uint i, j;
    for(i=pl->last - 1, j=0; j < pl->last; i=j++)
		v_sub(diffs + i, &pl->points[j].co, &pl->points[i].co);
    uint k, m;
    vec w;
    double d, ratio;
    for(i=pl->last - 1, j=0; j < pl->last; i=j++){
        for (k=(pl->last - 1), m=0; j < pl->last; k=m++){
            if (i == k || i == m || j == k || j == m)
                continue;
            d= v_perp(diffs + i, diffs + k, ax);
            if (fabs(d) < E)
                continue;
            v_sub(&w, &pl->points[i].co, &pl->points[k].co);
            ratio= v_perp(diffs + i, &w, ax) / d;
            if (0.0 < ratio && ratio < 1.0){
                ratio= v_perp(diffs + k, &w, ax) / d;
                if (0.0 < ratio && ratio < 1.0)
                    return 1;
            }
        }
    }
    return 0;
}

uint pl_simplify_2d(polygon pl, enum axis_pair ax)
{
    assert(pl);
    if (pl->last < 3)
        return pl->last;
    uint count= 0;
    uint i, j;
    vec diffs[pl->last];
    struct PL_Vertex *new_points= (struct PL_Vertex*)malloc(pl->last * sizeof(struct PL_Vertex));
    assert(new_points);
    for (i=(pl->last - 1), j=0; j < pl->last; i=j++)
        v_sub(diffs + i, &pl->points[i].co, &pl->points[j].co);
    for (i=(pl->last - 1), j= 0; j < pl->last; i=j++){
        if (fabs(v_perp(diffs + i, diffs + j, ax)) > ZEROAREA)
            new_points[count++]= pl->points[j];
    }
    free(pl->points);
    pl->points= new_points;
    pl->last= count;
    if (pl->holes){
        for (i=0; i < pl->holes->last; i++)
            pl_simplify_2d(pl->holes->polys[i], ax);
    }
    return count;
};

uint pl_simplify_3d(polygon pl)
{
    assert(pl);
    if (pl->last < 3)
        return pl->last;

	vec n;
    pl_normal(pl, &n);
	enum axis_pair ax;
	if (fabs(n.z) < ZEROLENGTH){
		if (fabs(n.y) > fabs(n.x)){
			ax= XZ;
		}
		else{
			ax= YZ;
		}
	}
	else{
		ax= XY;
	}
	return pl_simplify_2d(pl, ax);
};

uint pl_rm_doubles_2d(polygon pl, enum axis_pair ax)
{
    struct PL_Vertex *new_points= (struct PL_Vertex*)malloc(pl->last * sizeof(struct PL_Vertex));
    assert(new_points);
    uint i, j, count= 0;
    for (i= (pl->last - 1), j= 0; j < pl->last; i=j++){
        if (!same_2d(&pl->points[i].co, &pl->points[j].co, ax)){
            new_points[count++]= pl->points[j];
        }
    }
    if (count < pl->last){
        free(pl->points);
        pl->points= new_points;
        pl->last= count;
    }
    pl_recalc_bounds(pl);
    if (pl->holes){
        for (i=0; i < pl->holes->last; i++)
            pl_rm_doubles_2d(pl->holes->polys[i], ax);
    }
    return count;
}

uint pl_rm_doubles_3d(polygon pl)
{
    struct PL_Vertex *new_points= (struct PL_Vertex*)malloc(pl->last * sizeof(struct PL_Vertex));
    assert(new_points);
    uint i, j, count= 0;
    for (i= (pl->last - 1), j= 0; j < pl->last; i=j++){
        if (!same_3d(&pl->points[i].co, &pl->points[j].co)){
            new_points[count++]= pl->points[j];
        }
    }
    if (count < pl->last){
        free(pl->points);
        pl->points= new_points;
        pl->last= count;
    }
    pl_recalc_bounds(pl);
    if (pl->holes){
        for (i=0; i < pl->holes->last; i++)
            pl_rm_doubles_3d(pl->holes->polys[i]);
    }
    
    return count;
}

// keeps first point same
void pl_reverse(polygon pl)
{
    assert(pl);
    if (pl->last < 3)
        return; //nothing to reverse
    struct PL_Vertex *rpoints= (struct PL_Vertex*)calloc(pl->size, sizeof(struct PL_Vertex));
    struct PL_Vertex *last= pl->points + pl->last - 1;
    uint i;
    rpoints[0]= pl->points[0];
    rpoints[0].flags= last->flags;
    for (i=1; i < pl->last; i++){
        rpoints[i]= *last--;
        rpoints[i].flags= last->flags;
    }
    free(pl->points);
    pl->points= rpoints;
    if (pl->holes){
        for (i=0; i < pl->holes->last; i++)
            pl_reverse(pl->holes->polys[i]);
    }
}

void pl_prj_to_plane( polygon pl,
                      const vec * const pp, const vec * const pn,
                      const vec * const direction)
{
    assert(pl);
    assert(pn);
    assert(pp);
    assert(direction);
    double d= v_dot(pn, direction);
    assert(fabs(d) > E);// otherwise projection is impossible
    vec tmp, *pt;
    uint i;
    double N;
    for (i=0; i < pl->last; i++){
		pt= &pl->points[i].co;
        v_sub(&tmp, pt, pp);
        N= -v_dot(pn, &tmp);
        v_scale(&tmp, direction, N/d);
        v_add(pt, pt, &tmp);
    }
    v_sub(&tmp, &pl->bb.min, pp);
    N= -v_dot(pn, &tmp);
    v_scale(&tmp, direction, N/d);
    v_add(&pl->bb.min, &pl->bb.min, &tmp);
    
    v_sub(&tmp, &pl->bb.max, pp);
    N= -v_dot(pn, &tmp);
    v_scale(&tmp, direction, N/d);
    v_add(&pl->bb.max, &pl->bb.max, &tmp);
}

struct Dist {
    double distance;
    uint idx;
};

int cmp_dists(const void* dst_a, const void* dst_b)
{
    double d_a = ((struct Dist *)dst_a)->distance;
    double d_b = ((struct Dist *)dst_b)->distance;
    if (d_a < d_b)
        return -1;
    else if (d_a > d_b)
        return 1;
    else
        return 0;
}

void process_splits(polylist *splits_pt, enum PL_Order orig, enum axis_pair ax)
{
    assert(splits_pt);
    polylist splits= *splits_pt;
    assert(splits);
    polylist holes= pll_init(splits->last);
    polylist parts= pll_init(holes->size);
    polygon split;
    uint i;
    for (i=0; i < splits->last; i++){
        split= splits->polys[i];
        if (split->last < 3){
            pl_kill(&split);
            continue;
        }
        if (pl_order_2d(split, ax) == orig)
            pll_append(&parts, split);
        else
            pll_append(&holes, split);
    }
    

    if (holes->last){
        polylist leftovers, holes2join;
        uint j;
        polygon pl, hole;
        for (i=0; i < parts->last; i++){
            pl= parts->polys[i];
            holes2join= leftovers= NULL;
            for (j=0; j < holes->last; j++){
                hole= holes->polys[j];
                if (pl_rel_point_2d(pl, &hole->points[0].co, ax, 1) & OUT)
                    pll_append(&leftovers, hole);
                else if (hole->flag)
                    pll_append(&holes2join, hole);
                else
                    pl_add_hole(pl, hole);
            }

            if (holes2join){
                vec bpt= pl->points[pl->last - 1].co;
                uint limit= holes2join->last;
                struct Dist dists[limit];
                for (j=0; j < limit; j++){
                    hole= holes2join->polys[j];
                    dists[j].idx= j;
                    dists[j].distance= dist(&hole->points[0].co, &bpt);
                }
                qsort(&dists, limit, sizeof(struct Dist), cmp_dists);
                for (j=0; j < limit; j++){
                    hole= holes2join->polys[dists[j].idx];
                    pl_extend_verts(pl, hole->points, hole->last);
                }
                pll_killall(&holes2join); // points are safely transfered to part
            }
            free(holes->polys);
            free(holes);
            holes= leftovers;
            if (!holes)
                break; // all holes are distributed
        }
    }
    pll_kill(splits_pt);
    *splits_pt= parts;
}

int pl_split_by_plane(const polygon pl,
                      const vec * const pp, const vec * const pn,
                      polylist *fparts, polylist *bparts, int do_holes)
{
    assert(pl);
    uint limit= pl->last;
    uint i;
    int rmap[limit];
    int R= pl_rel_plane(pl, pp, pn, rmap);
    if (R==ON)
        return R;
    else if ((R == FR) || (R == ONFR)){
        if (fparts)
            pll_append(fparts, pl_copy(pl, NULL, 1));
    }
    else if((R == BK) || (R == ONBK)){
        if(bparts)
            pll_append(bparts, pl_copy(pl, NULL, 1));
    }
    else{

        int fcross, bcross;
        polylist last2first;
        polylist fsplits= pll_init(limit);
        polylist bsplits= pll_init(limit);

        bcross= ((*rmap == BK) || (*rmap == ON));
        fcross= !bcross;
        last2first= (fcross) ? fsplits : bsplits;

        uint j;
        int r, rprv, rnxt;
        vec isx;
        struct PL_Vertex * pt;
        polygon split= pl_init(limit);
        split->flag= 1;

        for (i=(limit - 1), j=0; j < limit; i=j++){
            pt= pl->points + i;
            r= rmap[i];
            rnxt= rmap[j];
            pl_append_vertex(&split, pt);
            if (r == ON){
                rprv= (i > 0) ? rmap[i-1] : rmap[limit - 1];
                R= rprv|rnxt;
                if (R == ON)
                    continue; //useless point ._*_.
                else if (R & BK){
                    if (++bcross % 2){
                        pll_append(&bsplits, split);
                        split= pl_init(limit);
                        split->flag= 1;
                        pl_append_vertex(&split, pt);
                        if (R == BK)
                            bcross++;
                    }
                }
                else if(R & FR){
                    if (++fcross % 2){
                        pll_append(&fsplits, split);
                        split= pl_init(limit);
                        split->flag= 1;
                        pl_append_vertex(&split, pt);
                        if (R == FR)
                            fcross++;
                    }
                }
            }
            if ((r|rnxt) == CROSS){
                line_intersect_plane(&pt->co, &pl->points[j].co,
                                     pp, pn, &isx);
                pl_append(&split, &isx, pt->flags);
                bcross++;
                fcross++;
                if (bcross % 2)
                    pll_append(&bsplits, split);
                else
                    pll_append(&fsplits, split);
                split= pl_init(limit);
                split->flag= 1;
                pl_append(&split, &isx, pt->flags);
            }
        }

        if (fcross % 2)
            pll_append(&fsplits, split);
        else
            pll_append(&bsplits, split);

        if (last2first->last > 1){
            polygon tojoin= last2first->polys[--last2first->last];
            pl_extend_verts(last2first->polys[0], tojoin->points, tojoin->last);
            pl_kill(&tojoin);
        }

        if (do_holes){
            if (pl->holes){
                for (i=0; i < pl->holes->last; i++)
                    pl_split_by_plane(pl->holes->polys[i],
                                      pp, pn, &fsplits, &bsplits, 0);
            }
            vec n;
            pl_normal(pl, &n);
            enum axis_pair ax;
            if (fabs(n.z) > ZEROLENGTH)
                ax= XY;
            else if (fabs(n.x) > fabs(n.y))
                ax= XZ;
            else
                ax= YZ;
            enum PL_Order ord= pl_order_2d(pl, ax);
            process_splits(&fsplits, ord, ax);
            process_splits(&bsplits, ord, ax);
        }
        // put parts back to the lists
        pll_extend(fparts, fsplits->polys, fsplits->last);
        pll_kill(&fsplits);
        pll_extend(bparts, bsplits->polys, bsplits->last);
        pll_kill(&bsplits);
    }
    return R;
}

//##################################### P O L Y L I S T ####
polylist pll_init(const uint count)
{
    polylist pll= (struct PL_List*)calloc(1, sizeof(struct PL_List));
    assert(pll);
    pll->polys= (struct Polygon**)malloc(count * sizeof(struct Polygon**));
    assert(pll->polys);
    pll->size= count;
    return pll;
}

polylist pll_check_pointer(polylist *pll_pt)
{
    assert(pll_pt);
    polylist result= *pll_pt;
    if (!result)
        result= *pll_pt= pll_init(10);
    return result;
}

void pll_expand(polylist pll, const uint count)
{
    assert(pll);
    pll->polys= (struct Polygon**)realloc(pll->polys, pll->size+= count);
    assert(pll->polys);
}

void pll_append(polylist *pll_pt, const polygon pl)
{

    polylist pll= pll_check_pointer(pll_pt);
    if(pll->last > pll->size){
        pll_expand(pll, 2 * pll->size);
    }
    pll->polys[pll->last++]= pl;
}

polygon pll_pop(polylist pll)
{
    assert(pll);
    return (pll->last) ? pll->polys[--pll->last] : NULL;
}

polylist pll_copy(polylist pll, void (*copy_method)(vec * dest, const vec * const src))
{
    if (!pll)
        return NULL;
    polylist cp= pll_init(pll->size);
    uint i;
    for (i=0; i < pll->last; i++)
        pll_append(&cp, pl_copy(pll->polys[i], copy_method, 1));
    return cp;
}

void pll_kill(polylist* pll_pt)
{
    assert(pll_pt);
    polylist pll= *pll_pt;
    assert(pll);
    if (pll->polys)
        free(pll->polys);
    free(pll);
    *pll_pt= NULL;
}

void pll_killall(polylist* pll_pt)
{
    assert(pll_pt);
    polylist pll= *pll_pt;
    assert(pll);
    if (pll->polys){
        uint i;
        for (i=0; i < pll->last; i++)
            pl_kill(&pll->polys[i]);
        free(pll->polys);
    }
    free(pll);
    *pll_pt= NULL;
}

int pll_remove_pl(polylist pll, const polygon pl)
{
    assert(pll);
    uint i;
    for (i=0; i < pll->last; i++){
        if (pll->polys[i] == pl)
            break;
    }
    if (i == pll->last)
        return 0;
    else
        pll_remove(pll, i);
    return 1;
}

polygon pll_remove(polylist pll, const uint idx)
{
    assert(pll);
    if (idx < pll->last){
        polygon result= pl_copy(pll->polys[idx], NULL, 1);
        uint i=idx;
        polygon* pl_pt= pll->polys + i++;
        for (; i < pll->last; i++)
            *pl_pt++= pll->polys[i];
        pll->last--;
        return result;
    }
    return NULL;
}

polygon pll_replace(polylist pll, const uint idx, const polygon pl)
{
    assert(pll);
    assert(pl);
    polygon result= pll->polys[idx];
    pll->polys[idx]= pl;
    return result;
}

int pll_replace_pl(polylist pll, const polygon to_replace, const polygon pl)
{
    assert(pll);
    uint i;
    for (i=0; i < pll->last; i++){
        if (pll->polys[i] == to_replace)
            pll_replace(pll, i, pl);
            return 1;
    }
    return 0;
}

void pll_extend(polylist *pll_pt, const polygon * const polys, const uint count)
{
    
    assert(polys);
    polylist pll= pll_check_pointer(pll_pt);
    if (!pll){
        *pll_pt= pll_init(2*count);
        pll= *pll_pt;
    }
    else if(pll->last + count > pll->size){
        pll_expand(pll, (pll->size * 2) + count);
    }
    uint i;
    for (i=0; i < count; i++)
        pll->polys[pll->last++]= polys[i];
}
 
uint pll_count_points(polylist pll)
{
    if(!pll)
        return 0;
    uint i, count=0;
    for (i=0; i < pll->last; i++)
        count+= pll->polys[i]->last;
    return count; 
}

int pl_cmp_area(const void *pl_pt_0, const void *pl_pt_1)
{
	double area_0= pl_area_3d(*((polygon *)pl_pt_0), 1);
	double area_1= pl_area_3d(*((polygon *)pl_pt_1), 1);
	if (area_0 < area_1)
		return 1;
	else if(area_0 > area_1)
		return 0;
	else
		return -1;
}
	
void pll_sort_by_area(polylist pll)
{
	qsort(pll->polys, pll->last, sizeof(polygon), pl_cmp_area);
}

